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ABSTRACT 

In this paper we present the results of numerical simulations of the Kelvin- 
Helmholtz instability in a stratified shear layer. This shear instability is believed to be 
responsible for extra mixing in differentially rotating stellar interiors and is the prime 
candidate to explain the abundance anomalies observed in many rotating stars. All 
mixing prescriptions currently in use are based on phenomenological and heuristic esti- 
mates whose validity is often unclear. Using three-dimensional numerical simulations, 
we study the mixing efficiency as a function of the Richardson number and compare 
our results with some semi-analytical formalisms of mixing. 
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1 INTRODUCTION 

Mixing is a fundamental process which profoundly affects 
the evolution of stars. Nonetheless, a theory of mixing that 
is applicable to stellar conditions is still missing. Hitherto 
stellar evolution theory has resorted to more or less heuris- 
tic prescriptions of mixing, the validity of which is often 
unclear. One of the prominent mechanisms responsible for 
mixing in stars is the shear or Kelvin-Helmholtz instability. 
This instability is particularly important since most stars 
are are known to be at least partly differentially rotating. 
Therefore, shear mixing is likely to play a significant role in 
stellar evolution. 

Observational evidence seems to suggest that current 
prescriptions underestimate the efficiency of the mixing pro- 
cesses at work, especially in fast rotating stars. The obser- 
vational evidence for mixing is increasing rapidly. Herrero et 
al. (1992) find that all fast rotating 0-stars show significant 
surface He-enrichments. This seems to suggest that mixing is 
strong enough to transport the nuclearly processed material 
to the surface in a fraction of the life time of the star. Other 
observations include the N/C and ^'^C and ^'^C enrichments 
of the Red Giant Branch (e.g. Kraft 1997, Charbonnel 1995), 
the He and N excesses in OBA supergiants (Fransson et al. 
1989), the depletion of boron in most B-type stars (Venn, 
Lambert & Lemke 1996) and the ratio of the number of 
blue to red supergiants in galaxies (Langer & Maeder 1995) . 
For a more detailed review of the observational evidence see 
Maeder (1995), Kraft (1994) and references therein. But the 
conclusion from these observations is that both, the enrich- 
ment in CNO elements and the depletion of fragile elements 
such as boron can be explained if some form of mixing is 



introduced (Langer 1997). 

However, a convincing theory of mixing under most as- 
trophysical conditions is still missing. One of the pioneering 
works in this field was performed by Endal & Sophia (1978) 
who estimated the efficiency of mixing due to several rota- 
tionally induced instabilities in massive stars. 

The Sun is also known to rotate differentially. Pinson- 
neault et al. (1989) treated the efficiency of rotational mix- 
ing as a free parameter and fitted it to a solar model. In the 
meantime, helioseismological observations from the space- 
craft SOHO and from the GONG network have yielded pre- 
cise measurements of the rotation rate in the solar interior. 
In a number of publications the effects of mixing on stellar 
evolution have been examined, see, e.g. Meynet & Maeder 
(1997) and Staritsin (1999). Recently, Denissenkov & Tout 
(2000) investigated the effects of 'deep mixing' on the evolu- 
tion of red giants. However, for want of a valid fundamental 
theory, the mixing formalisms employed in the aforemen- 
tioned studies have all been based on more or less heuristic 
estimates. It is the aim of this work to make steps towards 
a more fundamental theory of mixing. 

In chemically homogeneous, stratified shear flows, the 
Kelvin-Helmholtz instability occurs when the destabilising 
effect of the relative motion in the different layers domi- 
nates over the stabilising effect of buoyancy (see e.g. Chan- 
drasekhar 1961). The competition between the two effects 
is described by the Richardson number, Ri. For a non- 
dissipative, parallel, steady flow a simple linear stability cri- 
terion can be derived as follows: Let the horizontal velocity 
of the fluid be a function of z only, i.e. U = U{z). Now sup- 
pose that two volumes of fluid at heights z and z + Az are 
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interchanged where the densities at those two heights are p 
and p + Ap, and their horizontal velocities U and U + Af7, 
respectively. The work per unit mass that must be done 
against gravity is 



AVK = -gApAz. 



(1) 



After the exchange both cells have an average velocity of 
U = ^{2U + AU). Therefore, the kinetic energy available to 
do this work is 



\p[U'' + {U + AUf - i(2[/ + AUf] = y{AUf. 



Hence, a condition for stability is 



y{AUf < -gApAz, 



< 



^gdp 

p dz' 



(2) 



(3) 



(4) 



A more general criterion can be derived (see e.g. Shu 1992): 



Ri := 



> 



{dU/dzY ' 4' 
where N is the Brunt- Vaisala frequency given by 
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with 5 = —{din p/d\nT)p^^ and cj) = (91np/(?ln/u)p,T, 
where T is temperature, P is pressure and p, chemical po- 
tential. Hp is the pressure scale height, i.e. {din P/dz)^^ . 
Ri is the so-called Richardson number. Stability is predicted 
if the local Richardson number exceeds 1 /4 everywhere. 

More detailed analytical theories of shear mixing have 
been developed (Macdcr 1995, Canuto 1998) and some will 
be reviewed briefly in Sec. 2. In Sec. 3, we present the re- 
sults of three-dimensional simulations of shear instabilities 
in stratified fluids. These numerical simulations were per- 
formed for a simple configuration and a range of initial con- 
dition. Our aim is to elucidate the nonlinear dynamics of the 
shear instability and to study the parametric dependence of 
the mixing efficiency on the Richardson number. We quan- 
tify the mixing of the fluid by introducing 'tracer' particles 
that are advected with the fluid. With conditions in stellar 
interiors in mind, we only consider subsonic flows. From the 
motion of these particles a heuristic diffusion constant can 
be derived which is then compared to the analytic prescrip- 
tions described in Sec. 2. 

Apart from mixing by shear instabilities, there is an- 
other important candidate believed to be reponsible for ex- 
tra mixing in stars, namely convective overshoot. In mod- 
els of stellar evolution the fluid is assumed to bo mixed 
whereever some convective instability criterion (such as the 
Schwarzschild criterion) is fulfilled, and convectively sta^ 
ble where this criterion is not fulfilled. Mathematically the 




Figure 1 . Region of instability (shaded) in the K-Ri plane, where 
K is the dimensionless wave number and Ri the Richardson num- 
ber (see text). 



boundary between these two regions is sharp. In reality, how- 
ever, the convective fluid elements still carry some momen- 
tum when they reach this boundary and subsequently pene- 
trate into the convectively "stable" region. This mechanism 
is called convective overshoot and has been investigated nu- 
merically by a number of groups (Frcytag, Ludwig & Stef- 
fen 1996, Nordlund & Stein 1996, Singh, Roxburgh & Chan 
1998) all of which find some degree of overshoot. This con- 
clusion is supported by observations. The most compelling 
evidence stems from isochrone fitting to stellar clusters and 
binary systems which suggests that convective overshoot is 
significant (see Zahn 1991 for a review). 

Unlike convective overshoot which provides mixing at 
the boundaries of convection zones, we investigate mixing 
by rotational instabilities that operate in convectively stable 
regions. 



2 THEORIES OF SHEAR MIXING 

2.1 Linear stability analysis of a shear layer 

The stability of a stratified shear layer has boon studied ex- 
tensively using linear theory (sec e.g. Chandrasokhar 1961, 
Howard & Maslowe 1973, Maslowo & Thompson 1971, Hazel 
1972, Morris et al. 1990). Recently, Balmforth & Morri- 
son (1998) reviewed the conditions for instability in inviscid 
shear flows. In astrophysics, Kipponhahn & Thomas (1987) 
and MacDonald (1983) examined the redistribution of mat- 
ter and angular momentum in accreting white dwarfs by 
means of a linear stability analysis. In this section we will 
present a simple example of a linear stability analysis for a 
plane shear layer. The results of this analysis will then be 
compared against our numerical simulations in Sec. 3. 

The stability equation for wave-like perturbations in 
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its inviscid and incompressible form can be written as (e.g. 
Chandrasekhar 1961) 



U" 



U -c iU -cf 



= 0, 



(7) 



where (3 = p / p and cj>[z) &qi[ik{x — ct)\ is the perturba- 
tion stream function. Eqn. (^) is sometimes referred to as 
the Orr-Sommerfeld equation. As an example for the lin- 
ear analysis of the Kelvin-Helmholtz instability we will cite 
a case which was first considered by Taylor: It consits of 
two superposed fiuids of different densities separated by a 
transition layer of intermediate density in which the shear 
velocity varies linearly from that of the lower layer to that 
of the upper layer, i.e. 

layer 1: z > d, p = po(l — e), U — Uo 
layer 2: d > z > —d, p = po, U = Uoz/d 
layer 3: z < — d, p = po{l + e), U — —Uo- 

It is easily verified that the Richardson number for this 
problem is given by 



ged 



Ri=7^- (8) 

The solutions of eqn (^) in each of the three regions are 



layer 1 
layer 2 
layer 3 



(pi = Aie 
4>2 = A2e 
(f>3 = ylae* 



-kz 

-kz I D ^kz 
+ -fc>2e 



Continuity of (p at the interfaces, and requiring tp to 
vanish at ±oo, gives the characteristic equation: 



1 - 



Ri+ (£+ 1) + lek{c+l) 



1 - 



Ri- (c- 1) - iefc(c- 1)2 



(9) 



where k — 2kd and c — c/Uo- Eqn (9) can now be solved 
numerically for c for a range of wave and Richardson num- 
bers. When c becomes imaginary the perturbation becomes 
unstable. The parameter region in which this occurs, i.e. the 
region of instability, is shown in Fig. 1. 

The efficiency of mixing is sometimes parametrized by 
a diffusion coefficient which can be loosely defined as 



L»(Ri) 



Im 



.(Ri)] 



(10) 



where Im [cmax] is the maximal growth rate of the pertur- 
bations and fcmax the corresponding wave number - both 
evaluated at a given Richardson number. The diffusion co- 
efficient as a function of the Richardson number is shown 
in Fig. 2. As expected, D rapidly decreases with increasing 
Richardson number. 




Figure 2. Diffusion coefficient (in units of Uod) as a function of 
the Richardson number. 



(see e.g. Chandrasekhar 1961, Howard & Maslowe 1973) but 
qualitatively they all yield the same picture as shown in Fig. 
1. Linear stability analyses, however, are only valid for small 
perturbations. They can be useful for studying the onset of 
instability, but they are not very reliable in describing the 
form or the efficiency of mixing itself. This is a problem that 
requires numerical simulations and a set of such simulations 
will be presented below. 

In the next subsection we will quote a thermodynamic 
estimate of the diffusion coefficient, before we move on to 
the numerical simulations in Sec. 3. 



2.2 Maeder's prescription 

Maeder (1995, 1997) reexamined the Richardson criterion 
taking into account radiative losses (also see Maeder & Zahn 
1998). These authors derive a diffusion coefficient given by 



D = K 



{dU/dzf 
7V2 



(11) 



where K is the thermal diffusivity K = AacT^ /3Kp^Cp, a 
being the radiation density constant, c the speed of light, 
Cp the specific heat at constant pressure and k the opacity. 
Maeder's prescription for mixing due to differential rota- 
tion has been used, e.g., by Langer (1992) and Denissenkov 
(1994) for evolutionary calculations for O and B stars in 
our Galaxy and the Large Magellanic Cloud. More recently, 
Denissenkov & Tout (2000) invoked a diffusion constant of 
the form of Eq. (11) to explain extra mixing in globular 
cluster red giants. 



In this section we presented an example of just one of 
the many instability analyses found in the literature. The 
example was chosen for its simplicity. Of course, more realis- 
tic velocity and density distributions have been studied since 



3 NUMERICAL SIMULATIONS 

Here we present the results of numerical simulations of the 
dynamical shear instability in a stratified fluid. The simu- 
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Figure 3. Top: Initial density stratification; Middle: Shear ve- 
locities for some selected simulations; Bottom: Corresponding 
Richardson numbers. 

lations were obtained using the ZEUS-3D code which was 
developed especially for problems in astrophysical hydrody- 
namics (Clarke & Norman 1994). The code uses finite differ- 
encing on a Eulerian or pseudo-Lagrangian grid and is fiiUy 
explicit in time. It is based on an operator-split scheme with 
piecewise linear functions for the fundamental variables. The 
fluid is advcctcd through a mesh using the upwind, mono- 
tonic interpolation scheme of van Leer. For a detailed de- 
scription of the algorithms and their numerical implementar 
tion see Stone &i Norman (1992a, b). 

In our simulations we employed an ideal gas equation 
of state, we ignored the effects of magnetic fields, rotation, 
nuclear reactions and variations in radiative processes. The 
simulations were computed on a Cartesian grid and the com- 
putational domain was chosen to have the dimension 2 • 10*^ 
cm X 2 • 10* cm x 10* cm (in the x-, y- and z-direction, 
where gravity acts in the 2-direction) . It was covered by 100 
X 100 X 50 grid points. The calculations were performed on 
a CRAY Jedi parallel-processor and an IBM RS/6000 clus- 
ter. 

The gravitational acceleration was taken to be that of 
a point mass of 1 M© at a distance of 1 Rq from the 
lower boundary, where the gravitational acceleration acted 
purely in the ^-direction. An analytic model for an isother- 
mal density distribution under constant gravitational accel- 
eration was relaxed in ID until hydrostatic equilibrium was 




7.000>;10^° 7.010x10^'^ 7.000x10^° 7.010x10 

Figure 4. Radial and shear velocities at various times for the 
example with Ri=0.4 

attained. The residual velocities after the relaxation were 
six orders of magnitude smaller than the maximum shear 
velocities. The initial density distribution is shown in Fig. 
3. Then a shear velocity profile was imposed on the fluid. 
It was chosen to have the form of a hyperbolic tangent in 
order to minimise the effect of the boundaries onto the shear 
layer, i.e. 

U{z) = Uo tanh[(z - zo)/h], (12) 

where Uo is the amplitude of the shear velocity, z the vertical 
position of the shear layer, and h its extent. In order to keep 
the shear layer away from the boundaries, h was taken to 
be smaller than the vertical extent of the simulation region. 
Finally, Uo was chosen to yield a range of initial Richardson 
numbers of 0.05 - 3, where the Richardson number is taken 
in its original simple definition, i.e. Ri= gp' /pU'^, and is 
measured aX z = zq. Fig. 3 shows the initial conditions for 
Ri = 0.1, 0.4, 0.6 and 1.6. The boundary conditions were 
chosen to be periodic in the x and y direction and reflecting 
in the z direction. 

Figs. 4 and 5 show the x- and z-velocitics in a vertical 
slice through the computational domain at various times, 
for a flow with an initial Richardson number of 0.4 in the 
shear layer. The occurrence of turbulence is clearly seen. One 
can see that early on in the simulation the flow undulates 
slightly before the instability becomes fully nonlinear. Then 
one can observe that two vortices form which mix material 
over a large section of the computational domain. The vor- 
tices remain quite stable and increase in size. Their centres 
coincide with the centre of the shear layer. One can also 
note that at the centre of the vortex, the flow velocity is 
very small. For efficient transport the velocity fluctuations 
in the horizontal and the vertical have to be correlated. This 
correlation is apparent in the formation of vortices as seen 
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Velocity Field 




0.0 0,2 0.4 0.6 0.8 1.0 



-.=9000 s 



Velocity Field 




0.0 0.2 0.4 0.6 0.8 1.0 



-.= 12000 s 

Figure 5. Velocity field at t = 1500 s, t = 9000 s and t = 12000 
s for the same case as shown in FigJ4l 



in Fig. 5. More quantitatively it can be expressed through 
the vorticity. Maps of the vorticity in the direction both per- 
pendicular to the gravitational acceleration and the initial 
flow are shown in Fig. 6. (for Ri=0.4). The contours in Fig. 
6a. which show the vorticity after t — 1500 s span a range 
of 3 • 10~^ s~^ compared to a range of 10"^ s"^ in Fig. 6b. 
after 9000 s. This example shows the production of vorticity 
by the shear instability. 

The eddies have the form of long "rolls" which are 
aligned with the y-axis, i.e. perpendicular to, both, the shear 
flow and the gravitational acceleration. They are found to 
stay fairly symmetric in the j/-direction through the entire 
depth of the computational box. The structures on the big 
scales, as for example the "rolls" in our simulations, 'feel' the 
symmetry of the background model in the y-direction and 
therefore, they retain this symmetry. This has also been ob- 
served in laboratory experiments with incompressible fluids. 
Only on the small scales homogeneous, isotropic 3D tur- 
bulence will develop but we do not resolve these scales in 
our simulations. In the case of the shear instability, mixing 
is provided primarily by motions on big scales and, conse- 
quently, we argue that for our purpose we do not have to 
resolve these small scales. A further discussion of this issue 
can be found in Sec. 3.2. 

Furthermore, over longer distances the "rolls" will start 




Figure 6. Vorticity after t = 4500 s and t = 12000 s for the 
example shown in Figs. ^. In the top panel the contours span a 
range in vorticity of 3 ■ 10"'^ s~^ and in the bottom panel a range 
of 10-2 s-i. 




2x10^ 2x10^ 

Figure 7. Positions of the tracer particles for the example shown 
in Figs. 0. 



to feel the curvature of the star which is not accounted for 
in the slab geometry of our simulations and they may have 
the effect of breaking up the symmetry in the y-direction. 
These issues remain uncertainties in the results presented 
here and will have to be examined in a later study. 

In order to study mixing processes the ZEUS code was 
modified to follow the motion of 1000 'tracer' particles which 
are advected with the fluid. The positions of some of the 
tracer particles at various times (again for the case with Ri 
= 0.4) are shown in Fig. 7. 

The difi'usion constant can then be defined as 
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2.0x10^ 4.0x10^ e.OxlO'^ 8.0x10-^ 1.0x10* 1.2x10* 

lim=[s] 

Figure 8. Diffusion coefficient as a function of time for Ri=0.4. 
The solid line shows the result of the 3D simulation and the 
dashed line the result of the 2D simulation with the same resolu- 
tion. The dotted line represents the results from the computation 
on a 200 x 100 grid and the dash-dotted line on a 400 x 200 grid. 



D = a'^/t, (13) 

where '^'^ = ^^i\A^(^) ^ -^^o]^, A) being the original height 
of the A''th tracer particle and z{N) its height after a time t. 
As apparent from Fig. 7 the diffusion coefficient varies with 
time. The diffusion constant as a function of time have been 
plotted in Fig. 8. In Fig. 8 the Richardson number is 0.4, 
which is greater than 0.25 and, therefore, the flow should 
be stable to the shear instability according to the simple 
Richardson criterion. The simulations reveal that it is not, 
as has in fact been shown in previous simulations. One can 
observe that initially D rises with time before it eventually 
approaches a value which remains nearly constant for some 
time. The remaining scatter is mainly due to the stochastic 
nature of the turbulent mixing. The transient phase during 
which D rises is longer, the greater the Richardson number. 
Eventually, dissipation (mainly numerical) becomes notice- 
able and D starts to slowly decrease again. The constant 
value of D to which the curves are converging, is the value 
that is of greatest interest for the purpose of constructing 
stellar models. This value for D depends on the Richardson 
number. An understanding of the quantitative dependence 
of D on parameters such as the Richardson number will 
be useful for constructing stellar models which treat mixing 
through a diffusion equation. 

We have plotted this 'stationary' value of D as a func- 
tion of Ri for all the simulations which we have performed 
(see Fig. 10). The crrorbars indicate the residual scatter 
observed in the simulations. Configurations with Ri > 1.6 
proved to be stable for times up to 30,000 s. 

3.1 2D versus 3D simulations 

Two-dimensional simulations of convective instabilities are 
known to introduce numerical artefacts and to yield unre- 
alistic results. In 2D simulations of convection the energy 




0.0 1^ . . • -. . . . 1 
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t=8000 s 



Velocity Field 
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0.0 L ^ , . ^ . ^ . . , . . 

0.0 0.2 0.4 0,6 0.8 1.0 



t= 12000 s 

Figure 9. Velocity fields in a 2D simulation after t = s, t = 8000 
s and t = 12000 s. 



cascades to larger scales where viscous dissipation is less ef- 
fective, whereas 3D simulations reveal that the energy cas- 
cades to smaller scales where it is absorbed. Moreover, 2D 
simulations often produce 'rolls', i.e. vortices in the direction 
perpendicular to the simulated domain, which, for convec- 
tion, tend to overestimate the mixing. 

In order to investigate the importance of performing 

these simulations in three dimensions, the runs described 
above were repeated in two dimensions. The ZEUS code 
is designed to efficiently compute on 1- and 2-dimensionaI 
grids, so the runs could be readily repeated in 2D. The same 
initial model and the same set of physical parameters as in 
the 3D simulations were used, and a grid of 100 x 50 points 
was chosen. Fig. 9 shows a sequence of snapshots of the ve- 
locity. 
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0.2 0.4 0.6 0.8 
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Figure 10. Diffusion coefficients for different Richardson num- 
bers. 



The diffusion coefficient derived from the 2D simula- 
tions as a function of time is displayed as a dashed line in 
Fig. 8. The figures show that in 2D simulations D rises more 
slowly with time than in three dimensions. But eventually D 

converges to about the same value which was obtained iu the 
3D simulations. This is an encouraging result. Not only does 
this raise the trust in the credibility and self-consistency of 
these results, it also suggests that 2D simulations may be 
sufficient to determine the efficiency of shear mixing. This 
conclusion may not be too surprising given that the 3D sim- 
ulations revealed a flow structure consisting of long rolls that 
are parallel to the y-axis, i.e. perpendicular to the shear flow. 
On smaller scales, near the Kolmogoroff scale, the picture 
would look very different because one would observe the on- 
set of full isotropic 3D turbulence. But here we are only 
simulating the Kelvin-Helmholtz instability and not full 3D 
turbulence which is a much more difficult problem and which 
requires a much finer (and yet unattained) numerical reso- 
lution. This will be discussed further in the next section. 



3.2 Numerical viscosity 

One frequently voiced objection to these kinds of direct nu- 
merical simulations is that Reynolds numbers as high as 
those encountered in stars are unattainable on current com- 
puters, and that therefore the results are unrealistic. But, 
as pointed out by Balbus, Hawley & Stone (1996), this crit- 
icism is unjustified for simulations of the shear instability. 
They argue, that in order to simulate the onset of instability 
in a laminar flow, it is merely necessary that the 'typical' 
wavelength of the instability is resolved by the numerical 
scheme and that the numerical diffusion at this wavelength 
is less than the growth rate. This makes the simulation of 
shear instabilities an easier task than the simulation of vis- 
cous instabilities where in theory one would have to resolve 



everything down to the viscous length scale. 

If the flow is unbounded, there is no viscous boundary 
layer which might interfere with the results. The nonlinear 
instabilities are fundamentally inviscid in character. There- 
fore, we only need a resolution capable of resolving a range 
of wavelengths with numerical diffusion errors less than the 
growth rates. This view is also supported by experimental 
observations which suggest that mixing through shear instar 
bilities is dominated by large scale motions. 

In Fig. 12 we have plotted the kinetic energy in the 

vertical (^-direction) as a function of time. The vertical ki- 
netic energy is expressed in units of the total initial kinetic 
energy. In the beginning the energy is very small until the 
instability sets in. Then the kinetic energy jumps to about 
10% of its total initial value and slowly decreases amidst 
big fluctuations. This decrease is mainly due to numerical 
viscosity and shows yet again that we are not simulating 
proper 3D turbulence which would show a saturation of the 
kinetic energy. 

Porter & Woodward (1994) have estimated the 
Reynolds number of hydrodynamical simulations based on 
a PPM (piecewise parabolic method) code. The Reynolds 
number depends on the truncation error of the finite- 
difference algorithm, the Courant number, the background 
adveetion and the rmmber of grid points. They found that 
the effective Reynolds number is proportional to the third 
power of the number of grid points, with the main dissi- 
pation occurring at short wavelengths. Above this critical 
wavelength the diffusion was found to be small. Since the 
ZEUS code uses piecewise linear functions instead of piece- 
wise parabolic functions, the truncation errors of ZEUS will 
be larger than those of a PPM code. But even if they were 
ordy proportional to the second power of the number of grid 
points in one direction, we would still expect a rmmerical 
Reynolds number of around 10**. 

The effect of the numerical viscosity is to suppress fluc- 
tuations beneath the 'effective' viscous length scale. This 
will become important for high enough Richardson num- 
bers and all numerical simulations will find stability above 
a certain Richardson number. In our ease the fiow remained 
stable for Ri >1.6. Higher resolutions, and therefore higher 
Reynolds numbers will push this limiting Richardson num- 
ber to slightly higher values. 

3.3 Dependence on grid size 

We have investigated the dependence of our results on the 
resolution of the computational grid by repeating the 2D 
runs presented above (100 x 50) on grids with (b) 200 x 
100 and (c) 400 x 200 grid points. The variation of D with 
time is shown in Fig. 8. by the dotted (b) and dash-dotted 
(c) lines. It is immediately evident that the behaviour of 
D with time varies with the numerical resolution. For more 
finely resolved grids, D rises more quickly. This was to be 
expected since higher resolutions produce higher Reynolds 
numbers and therefore a less viscous flow. For higher effec- 
tive Reynolds numbers the turbulence grows more rapidly, 
and this is observed here. But again it merits mention that 
for all simulations D eventually converges to rougly the same 
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value for a given Richardson number, as long as the Richard- 
son number is not too big. For large Richardson numbers, 
the numerical viscosity becomes important, but here we will 
restrict the discussion to Richardson numbers less than 2. 



Finally, it can be seen in Fig. 5 that the eddies nearly 
span the entire vertical extent of the simulation domain. 
Hence, it is interesting to verify whether the vertical size 
of the computational domain affects the size of the eddies, ^ 

and therewith the value of the mixing coefficients. Wc have ^ 
repeated some 2D simulations in computational domains ^ 
which had a vertical extent of more than 50h or ten times "n" 
the vertical size of the previously shown simulations. Wc did 
not find that the eddies became bigger when the size of the 
computational box was increased. 



4 DISCUSSION 

In this paper results of direct numerical simulations of shear 
instabilities were presented. The simulations were carried 
out using the ZEUS hydro code. A chemically homogeneous 
stratified fluid in hydrostatic equilibrium was set up, onto 
which a shear flow was superimposed. The mixing was quan- 
tified by studying the dispersion of 'tracer' particles that are 
advcctcd by the fluid. Thus, the variation of the diffusion 
constant with the Richardson number was examined. 

In accordance with previous simulations, it was found 
that efficient mixing occurs for Richardson numbers sub- 
stantially higher than 0.25, contrary to the simple linear 
stability criterion. Now the diffusion coefficients that were 
derived from the motion of the tracer particles can be com- 
pared to the analytical estimates discussed in Sec. 2. Apply- 
ing Maeder's formalism to our initial model, i.e. Eqn (11), 
yields substantially lower values for D than found in our 
simulations. For our initial models, Maeder's diffusion coef- 
ficients in the centre of the shear layer lie between values of 
10® - 10^ cm /s. This is more than 3 orders of magnitude 
smaller than the diffusion coefficients found here. Finally, 
one might want to compare our results with Fig. 2, even 
though this figure assumed a different shear velocity distri- 
bution than the hyperbolic tangens used in the simulations. 
In Fig. 11 where the dimensionless diffusion coefficient (in 
units of Uoh) is plotted against the Richardson number. This 
dimensionless presentation is also useful for a wider appli- 
cation of our results. It is interesting to note that now our 
values for D are smaller by almost two orders of magnitude 
than those predicted by the linear analysis of our simple 
model of Sec. 2. 

It was argued that numerical viscosity has no great ef- 
fect on our findings Eis the effective Reynolds number is small 
on the scales on which the Kolvin-Helmholtz instability op- 
erates. This is true as long as the Richardson number is 
not too big. Furthermore, the differences between 2D and 
3D simulations of shear instabilities were studied. We found 
that the dimensionality of the simulation hardly affects the 
final 'steady state' diffusion coefficient. It only affects the dy- 
namics of the onset of the instability until a quasi-stationary 
state has been reached. The same qualitative observation 
was made when the resolution of the numerical grid was var- 
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Figure 11. Diffusion coefficients for different Richardson num- 
bers in units of Uoh (errorbars have been omitted). 

ied. It was found that the resolution affects the initial tem- 
poral variation of the diffusion coefficient, but hardly affects 
the final value of D. However, the similarity between the 2D 
and 3D simulations is only valid on the big scales which we 
arc resolving. On smaller scales, near the Kolmogoroff scale, 
one would observe full isotropic 3D turbulence which looks 
very different in two dimensions. 

To make the Richardson number the only controlled pa- 
rameter in our study, is clearly a simplification of the factors 
that determine the efficiency of mixing. In reality, the effi- 
ciency of mixing will depend on the density stratification 
and the velocity gradient separately, and not solely on the 
Richardson number. Moreover, the mixing will depend on 
the exact shape of the velocity profile and, only to a first 
approximation, on its first derivative. Nevertheless, if one is 
interested in the efficiency of mixing as a function of a single 
parameter, e.g. for stellar evolution studies, the Richardson 
number is still the most suitable parameter, since in the limit 
of an inviscid and incompressible shear flow, the stability is 
governed by the Richardson number. 

In this work the effects of the dynamical shear in- 
stability in a plane shear layer were investigated. Clearly, 
this is not the only mechanism responsible for mixing in 
stars. In rotating fluid bodies other instabilities may occur, 
such as, e.g., the Solberg-Holand instability, the Eddington- 
Swcct instability and the Goldroich-Schubcrt-Fricko instabil- 
ity. However, the shear instability operates on the dynamical 
timescale whereas the instabilities tied to the rotation occur 
on the longer Eddington-Vogt timescale. Therefore, these in- 
stabilities operate on a much longer timescale and therefore 
play a lesser role in diffusion processes, even though they 
can affect the long-term evolution of the star. 

The consequences of the newly found diffusion coeffi- 
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time 

Figure 12. Kinetic energy in the vertical z-direction as a function 
of time. The kinetic energy is expressed in units of the total initial 
kinetic energy and the units of time are arbitrary. 

cients onto stellar evolution and surface abundances are dif- 
ficult to foresee. In stars the speed as well as the depth of 
mixing determine the balance between mixing and nuclear 
burning. Therefore, the effect of mixing depends sensitively 
on the conditions prevailing in the star, as for example the 
positions and extents of its convection and burning shells. 

Finally we should mention that a number of factors can 
inhibit or facilitate mixing such as gradients in the chemical 
potential of the fluid, diffusion of radiation, magnetic fields 
and effects pertaining to the spherical geometry. In late-type 
stars strong chemical composition gradients exist, which will 
have a stabilising effect on the stratification. Therefore, espe- 
cially for stellar evolution studies, the chemical composition 
gradient is an important parameter, which will have to be 
included in future work. Similarly, poloidal magnetic fields, 
as they are believed to exist in the Sun, will also suppress 
mixing processes. These effects are potentially of great im- 
portance in stars and will be studied in a forthcoming paper. 
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